MODULE mod_chem_spack_qssa

  USE mem_chem1, ONLY: &
       chem1_vars          ! Type

  USE spack_utils, ONLY: &
       index_dim           ! Type

  USE rconstants, ONLY: &
       cp,              & ! PARAMETER
       cpor,            & ! PARAMETER
       p00                ! PARAMETER

  USE mod_chem_spack_kinetic, ONLY: &
       kinetic  ! Subroutine
  
  USE mod_chem_spack_fexloss, ONLY: &
       fexloss  ! Subroutine

  USE mod_chem_spack_fexprod, ONLY: &
       fexprod  ! Subroutine

  USE mod_chem_spack_dratedc, ONLY: &
       dratedc  ! Subroutine

  USE mod_chem_spack_rates, ONLY: &
       rates    ! Subroutine

  USE chem1_list, ONLY: nr

  USE FastJX, ONLY: &
       fast_JX_vars      ! Type
  
  USE FastJX, ONLY:  &
       fast_JX_g        ! %jphoto(IN)

   USE uv_atten, ONLY: &
       uv_atten_g 

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: chem_qssa ! Subroutine
 

CONTAINS


  !========================================================================================
  SUBROUTINE chem_qssa(m1,m2,m3,itchim,nrmax,mynum,ngrid,dtlt,pp,pi0,theta,rv,maxgrds, &
                       cosz,fjphoto,nspecies,nr_photo,weight,PhotoJMethod,chem1_g,index_g, &
                       nob,maxblock_size)
  !========================================================================================

    ! mem_grid
    INTEGER          , INTENT(IN)    :: ngrid
    REAL             , INTENT(IN)    :: dtlt

    ! node_mod
    INTEGER          , INTENT(IN)    :: mynum
    INTEGER          , INTENT(IN)    :: m1
    INTEGER          , INTENT(IN)    :: m2
    INTEGER          , INTENT(IN)    :: m3
    INTEGER          , INTENT(IN)    :: maxgrds

    ! mem_basic
    REAL             , INTENT(IN)    :: pp(m1,m2,m3)
    REAL             , INTENT(IN)    :: pi0(m1,m2,m3)
    REAL             , INTENT(IN)    :: theta(m1,m2,m3)
    REAL             , INTENT(IN)    :: rv(m1,m2,m3)

    ! chemistry_driver() local variable
    INTEGER          , INTENT(IN)    :: itchim(maxgrds)

    ! mem_radiate
    REAL             , INTENT(IN)    :: cosz(m2,m3)

    ! chem1_list
    INTEGER          , INTENT(IN)    :: nspecies
    INTEGER          , INTENT(IN)    :: nr_photo
    INTEGER          , INTENT(IN)    :: nrmax
    REAL             , INTENT(IN)    :: weight(nspecies)
    CHARACTER(LEN=10), INTENT(IN)    :: PhotojMethod

    ! FastJX
    TYPE(fast_JX_vars), INTENT(IN)    :: fjphoto

    ! mem_chem1
    TYPE (chem1_vars), INTENT(INOUT) :: chem1_g(nspecies,maxgrds)

    ! spack_utils
    TYPE(index_dim) ,   INTENT(IN)    :: index_g(maxgrds)
    INTEGER         ,   INTENT(IN)    :: nob(maxgrds)
    INTEGER         ,   INTENT(IN)    :: maxblock_size


    !  INTEGER,PARAMETER :: block_size=45
    !  INTEGER,PARAMETER :: block_size=4
    !  INTEGER,PARAMETER :: block_size=1
    DOUBLE PRECISION,PARAMETER     :: pmar=28.96  
    ! !LPCE
    !srf
    INTEGER :: i,ijk,n,cont
    DOUBLE PRECISION dtlt_i,fxc
    INTEGER k_,i_,j_,kij_
	INTEGER :: ii,jj


    TYPE spack_type_
       !3d Real
       DOUBLE PRECISION,POINTER,DIMENSION(:,:,:) :: dw
       !2d Real
       DOUBLE PRECISION,POINTER,DIMENSION(:,:)   :: jphoto 
       DOUBLE PRECISION,POINTER,DIMENSION(:,:)   :: loss	
       DOUBLE PRECISION,POINTER,DIMENSION(:,:)   :: prod	
       DOUBLE PRECISION,POINTER,DIMENSION(:,:)   :: rk    
       DOUBLE PRECISION,POINTER,DIMENSION(:,:)   :: w    
       DOUBLE PRECISION,POINTER,DIMENSION(:,:)   :: sc_p	
       !1D Real
       DOUBLE PRECISION,POINTER,DIMENSION(:)     :: temp	
       DOUBLE PRECISION,POINTER,DIMENSION(:)     :: press 
       DOUBLE PRECISION,POINTER,DIMENSION(:)     :: cosz	  
       DOUBLE PRECISION,POINTER,DIMENSION(:)     :: att	  
       DOUBLE PRECISION,POINTER,DIMENSION(:)     :: vapp 
       DOUBLE PRECISION,POINTER,DIMENSION(:)     :: volmol 
       !DOUBLE PRECISION,POINTER,DIMENSION(:)     :: xlw 
    END TYPE spack_type_

    TYPE(spack_type_) :: spack

!!$  !Determining number of blocks and masks
!!$  CALL AllocIndex(block_size,mynum,mmxp,mmyp,mmzp,mia,miz,mja,mjz,mibcon,mi0,mj0,&
!!$       ngrids,dtlongn)

    dtlt_i=1./dtlt

    !- Allocating Spaces to copy structure
    !ALLOCATE(spack(nob))
    !- 3D Variables
    ALLOCATE(spack%dw    (1:maxblock_size,nrmax,nspecies))
    !- 2D Variables
    ALLOCATE(spack%jphoto(1:maxblock_size,nr_photo))	;spack%jphoto=0.d0
    ALLOCATE(spack%sc_p  (1:maxblock_size,nspecies))	;spack%sc_p=0.d0
    ALLOCATE(spack%loss  (1:maxblock_size,nspecies))	;spack%loss =0.d0
    ALLOCATE(spack%prod  (1:maxblock_size,nspecies))	;spack%prod  =0.d0
    ALLOCATE(spack%rk    (1:maxblock_size,nrmax))		;spack%rk   =0.d0
    ALLOCATE(spack%w     (1:maxblock_size,nrmax))		;spack%w    =0.d0
    !- 1D variables
    ALLOCATE(spack%temp  (1:maxblock_size))			;spack%temp =0.d0
    ALLOCATE(spack%press (1:maxblock_size))			;spack%press =0.d0
    ALLOCATE(spack%cosz  (1:maxblock_size))			;spack%cosz  =0.d0
    ALLOCATE(spack%att   (1:maxblock_size))			;spack%att   =0.d0
    ALLOCATE(spack%vapp  (1:maxblock_size))			;spack%vapp =0.d0
    ALLOCATE(spack%volmol(1:maxblock_size))			;spack%volmol=0.d0
    !ALLOCATE(spack%xlw   (1:maxblock_size))

    DO i=1,nob(ngrid)

       !- Copying structure from input to internal
       DO ijk=1,index_g(ngrid)%block_end(i)

          k_=index_g(ngrid)%indexk(ijk,i) ! k brams index
          i_=index_g(ngrid)%indexi(ijk,i) ! i brams index
          j_=index_g(ngrid)%indexj(ijk,i) ! j brams index

          !- Exner function divided by Cp
          spack%press(ijk)=(  pp (k_,i_,j_) &
               + pi0(k_,i_,j_) )/cp
          !- Air temperature (K)			      			      
          spack%temp(ijk)= theta(k_,i_,j_)*spack%press(ijk)

          !- transform from Exner function to pressure
          spack%press(ijk)= spack%press(ijk)**cpor*p00

          !- Water vapor 
          spack%vapp(ijk)=  rv(k_,i_,j_)!&
          !* basic_g(ngrid)%dn0(k_,i_,j_)*6.02e20/18.


          !spack%volmol(ijk)=(spack%press(ijk)*1.e-3)/(8.314e0*spack%temp(ijk))
          spack%volmol(ijk)=(spack%press(ijk))/(8.314e0*spack%temp(ijk))
          !LPCE
          DO n=1,nspecies

             !LPCE conversion from ppbm to molecule/cm^3
             !spack%sc_p(ijk,n)=chem1_g(n,ngrid)%sc_p(k_,i_,j_) &
             !                    *(6.02e23*1e-12*spack%volmol(ijk))
             spack%sc_p(ijk,n)=chem1_g(n,ngrid)%sc_p(k_,i_,j_) &
                  *(6.02e23*1e-15*spack%volmol(ijk))*pmar/weight(n)
!	print *, chem1_g(n,ngrid)%sc_p(k_,i_,j_),spack%sc_p(ijk,n);call flush(6)
          END DO
          !spack%xlw(ijk) =  micro_g(ngrid)%rcp(k_,i_,j_) &
          !                   * basic_g(ngrid)%dn0(k_,i_,j_)*1e-3
          spack%cosz(ijk)=cosz(i_,j_)
       END DO

       !- FAST-JX section
    IF(trim(PhotojMethod) == 'FAST-JX' .or. trim(PhotojMethod) == 'FAST-TUV' ) THEN
          DO ijk=1,index_g(ngrid)%block_end(i)

             k_=index_g(ngrid)%indexk(ijk,i) ! k brams index
             i_=index_g(ngrid)%indexi(ijk,i) ! i brams index
             j_=index_g(ngrid)%indexj(ijk,i) ! j brams index

             DO n=1,nr_photo
                spack%jphoto(ijk,n)= fjphoto%jphoto(n,k_,i_,j_)
             END DO

          ENDDO
     
       ELSEIF(trim(PhotojMethod) == 'LUT') then
      
          DO ijk=1,index_g(ngrid)%block_end(i)
        
            i_=index_g(ngrid)%indexi(ijk,i) ! i brams index
            j_=index_g(ngrid)%indexj(ijk,i) ! j brams index
        
	   !- UV  attenuation un function AOT
           spack%att(ijk)= uv_atten_g(ngrid)%att(i_,j_) 
        
	   !- get zenital angle (for LUT PhotojMethod)
            spack%cosz(ijk)=  cosz(i_,j_)

          ENDDO     
       ENDIF

       !- call kinetic
       ! IN: nr_photo,spack%Jphoto,spack%temp,spack%vapp,spack%Press,spack%cosz
       !     1.0D0,1,index_g(ngrid)%block_end(i),maxblock_size,nrmax
       ! OUT: spack%rk
       CALL kinetic(nr_photo,spack%Jphoto   &
                           ,spack%rk    &
                           ,spack%temp  &
			   ,spack%vapp   &
                           ,spack%Press &
		           ,spack%cosz  &
    			   ,spack%att     &
			   ,1,index_g(ngrid)%block_end(i),maxblock_size,nr)
 
       CALL aqss(DTLT,spack%RK,spack%sc_p,spack%w,spack%dw,nr,nspecies,1,&
                 index_g(ngrid)%block_end(i),maxblock_size)
!print *,"chamou aqss",i
       !- restoring species loss/prod from internal to ccatt-brams structure
       DO ijk=1,index_g(ngrid)%block_end(i)
          !- tendency index (mude para calcular somente uma vez, isto e', crie um kij(ijk,i)  )
          k_   =index_g(ngrid)%indexk(ijk,i)   !- k   brams index
          i_   =index_g(ngrid)%indexi(ijk,i)   !- i   brams index
          j_   =index_g(ngrid)%indexj(ijk,i)   !- j   brams index
          kij_   =index_g(ngrid)%kij_index(ijk,i)!- kij brams tendency index

          DO n=1,nspecies
             !LPCE conversion from molecule/cm^3 to ppbm
             !chem1_g(n,ngrid)%sc_p(k_,i_,j_)= spack%sc_p(ijk,n)&
             !                  /(6.02e23*1e-12*spack%volmol(ijk))
             fxc = 1./(6.02e23*1e-15*spack%volmol(ijk)*(pmar/weight(n)))
             chem1_g(n,ngrid)%sc_t(kij_) = chem1_g(n,ngrid)%sc_t(kij_)   +  &    ! previous tendency
                  (  spack%sc_p(ijk,n) * fxc  -  &    ! new mixing ratio
                  chem1_g(n,ngrid)%sc_p(k_,i_,j_))&! old mixing ratio
                  *  dtlt_i                              ! inverse of timestep

          END DO
       END DO


    END DO

    !Deallocating local structure for later use
    !3D Variables
    DEALLOCATE(spack%dw)
    !2D Variables
    DEALLOCATE(spack%jphoto)
    DEALLOCATE(spack%sc_p)
    DEALLOCATE(spack%loss)
    DEALLOCATE(spack%prod)
    DEALLOCATE(spack%rk)
    DEALLOCATE(spack%w)
    DEALLOCATE(spack%cosz)
    DEALLOCATE(spack%att)
    !1D variables
    DEALLOCATE(spack%temp)
    DEALLOCATE(spack%press)
    DEALLOCATE(spack%vapp)
    DEALLOCATE(spack%volmol)
    !DEALLOCATE(spack%xlw)
    !DEALLOCATE(spack)

!!$  !Deallocating blocks data
!!$  CALL AllocIndex(block_size,mynum,mmxp,mmyp,mmzp,mia,miz,mja,mjz,mibcon,mi0,mj0,&
!!$         ngrids,dtlongn)

  END SUBROUTINE chem_qssa

SUBROUTINE aqss(DLTT,RK,Y,W,DW,nr,nspecies,ijkbegin,ijkend,maxblock_size)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!!
!	Avanco temporal em DLTT do sistema de EDO's Stiff pelo metodo alfa-QSS (mott99)
!	
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! eps  : convergencia geral
! eps1 : escala do passo temporal
! itmax :numero maximo de iteracoes do corretor
! iter : iteracao corrente
! DLTT : passo global de integracao
! ijkbegin e ijkend : inicio e fim de bloco (posic espacial)
! nr   : numero de reacoes do modelo cinetico
! ngas : numero de especies quimicas consideradas no modelo cinetico
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!	Declaracao de variaveis
!
implicit none
INTEGER, INTENT(IN) :: nspecies,nr,ijkbegin,ijkend,maxblock_size
INTEGER :: i,ijk,ITER, j,CONT,ITER_C,CONTMAX
REAL :: DLTT
DOUBLE PRECISION, external :: maximo,minimo,Pade

DOUBLE PRECISION :: DT,EPS2,RSIG,SIGMA,T,x
DOUBLE PRECISION :: Ymin,F,FC
DOUBLE PRECISION,ALLOCATABLE :: A(:,:),DTcel(:),E(:),G0(:,:),P0(:,:),P(:,:),Q(:,:),Q0(:,:),R(:), &
                                SIGMAcel(:),Y0(:,:),Yp(:,:)
DOUBLE PRECISION, INTENT(IN) :: RK(:,:)
DOUBLE PRECISION, INTENT(INOUT) :: Y(:,:),DW(:,:,:),W(:,:)
! DOUBLE PRECISION,ALLOCATABLE :: 	
integer, parameter :: ITMAX = 3
real, parameter :: eps = 0.01, eps1 = .01

ALLOCATE (DTcel(1:maxblock_size))
ALLOCATE (SIGMAcel(1:maxblock_size))
ALLOCATE (A(1:maxblock_size,1:nspecies))
ALLOCATE (E(1:nspecies))
ALLOCATE (G0(1:maxblock_size,1:nspecies))
ALLOCATE (P0(1:maxblock_size,1:nspecies))
ALLOCATE (P(1:maxblock_size,1:nspecies))
ALLOCATE (Q(1:maxblock_size,1:nspecies))
ALLOCATE (Q0(1:maxblock_size,1:nspecies))
ALLOCATE (R(1:nspecies))
ALLOCATE (Y0(1:maxblock_size,1:nspecies))
ALLOCATE (Yp(1:maxblock_size,1:nspecies))

!DLTT = 600.

Ymin = 0.
F = 10.				!Parametro que decide se Qo << PYo na inicializacao
EPS2 = EPS/1.5								! eq. (3.36)
!
!	Taxas de reacoes elementares e derivada
!
	CALL cinetica(RK,Y,W,DW,P,Q,nspecies,ijkbegin,ijkend,maxblock_size,nr)
!
!	Inicializacao e parametros do esquema
!
T = 0.
CONT = 0
CONTMAX = 5
ITER = 0
ITER_C = 0
!			print *,"************************************************************************************"
!			print *, " 					Composicao inicial"
!		do ijk = 1, maxblock_size
!			print *,"************************************************************************************"
!			print *, " "
!			print *, "bloco ", ijk
!			print *, " "
!			print *, "***********************************************************************************"
!			do j=1,nspecies
!				print*,"especie",j,"Y=", Y(ijk,j)				
!			enddo
!		enddo

do while(T.lt.DLTT)
!do while(CONT.LE.1000)
	if(ITER.eq.0) then	
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%						
!
!	 								INICIALIZACAO		
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%				
		do ijk=ijkbegin,ijkend	
!
! 			Calculo de G0 e R													(eq. 3.40)
!
			do j = 1,nspecies
				Y0(ijk,j) = Y(ijk,j)			! Salva uma copia de Y0
				P0(ijk,j) = P(ijk,j)			! Salva uma copia de P0	
				Q0(ijk,j) = Q(ijk,j)			! Salva uma copia de Q0							
				G0(ijk,j) = Q(ijk,j) -  P(ijk,j)*Y(ijk,j) 
				if (Q(ijk,j).gt.(F*P(ijk,j)*Y0(ijk,j)).and.P(ijk,j).ne.0.) then	
					R(j) = 1./P(ijk,j)
				elseif (Y(ijk,j).ne.0.) then
					R(j) = Y(ijk,j)/(dabs(G0(ijk,j)) + 1.d-6) ! Caso a taxa dy/dt seja
									      ! zero (equilibrio) penalizar Rj
				else
					R(j) = 1.d6				  
   				endif
!				print *, j,Y(ijk,j),Q(ijk,j),P(ijk,j),G0(ijk,j),R(j); call flush(6)
			end do
!
!	Passo temporal inicial
!
			DTcel(ijk) = eps1*minval(R)
!			DTcel(ijk) = eps1*minimo(R,nspecies)
!			print *, ijk,"calculou DT inicial", DTcel(ijk); call flush(6)	
		enddo

	DT = minval(DTcel)
!	DT = minimo(DTcel,maxblock_size)
	ITER = ITER + 1
	endif
!	print * ,"passo temporal do preditor", DT,"DLTT =", DLTT
!	print *, " "
!	print *, "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"
!	print *, " "
!	call flush(6)
				
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! 								P R E D I T O R 															
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if(ITER_C.eq.0) then
		do ijk=ijkbegin,ijkend
			do j=1,nspecies
				if(P0(ijk,j).ne.0.) then
				       x=1./(P0(ijk,j)*DT)
				       A(ijk,j)=	 (180.*x*x*x + 60.*x*x + 11.*x + 1.)/(360.*x*x*x + 60.*x*x + 12.*x + 1.) 

					!A(ijk,j) = Pade(1./(P0(ijk,j)*DT))				! A0
				else
					A(ijk,j) = .5
				endif
				Yp(ijk,j) = Y0(ijk,j)+DT*G0(ijk,j)/(1.+A(ijk,j)*DT*P0(ijk,j))				
			end do
		enddo												! ijk > bloco
	endif													! se corretor não iterou	
	
	CALL cinetica(RK,Yp,W,DW,P,Q,nspecies,ijkbegin,ijkend,maxblock_size,nr)
	
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!								Interpolacoes
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
!	Calculo de P barra, A barra e Q til	usando Yp													(eq. 3.22) 
!
	do ijk=ijkbegin,ijkend		
		do j=1,nspecies
			P(ijk,j) = .5*(P0(ijk,j) + P(ijk,j))					! P barra
			if(P(ijk,j).ne.0.) then
			      x=1./(P0(ijk,j)*DT)
			      A(ijk,j)= 	(180.*x*x*x + 60.*x*x + 11.*x + 1.)/(360.*x*x*x + 60.*x*x + 12.*x + 1.) 
                             !A(ijk,j) = Pade(1./(P0(ijk,j)*DT))			       ! A0
			else
				A(ijk,j) = .5
			endif

			Q(ijk,j) = A(ijk,j)*Q(ijk,j) + (1. - A(ijk,j))*Q0(ijk,j)		! Q agora contem Q til			(eq. 3.32)
		end do
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! 								C O R R E T O R 															
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		do j=1 ,nspecies
			Y(ijk,j) = Y0(ijk,j)+DT*(Q(ijk,j)-P(ijk,j)*Y0(ijk,j))/(1.+A(ijk,j)*DT*P(ijk,j))	! Y contem Yc		(eq. 3.31)
!	print *, ijk,j,Y0(ijk,j),Yp(ijk,j),Y(ijk,j); call flush(6)
!
!	Verificar quais especies atingiram a concentracao minima
!
! 	Variacao relativa de Y para o corretor											(eq. 3.34)
!
   			if (Y(ijk,j).gt.Ymin) then 					! verificar se Ymin = threshold da chem_spack_gear
   				E(j) = abs(Y(ijk,j) - Yp(ijk,j))/Y(ijk,j)		
			else
				E(j) = 0.
   			endif
		end do
		
		SIGMAcel(ijk) = maxval(E)/EPS2
!		SIGMAcel(ijk) = maximo(E,nspecies)/EPS2
!		print *, "celula",ijk,"sigma",SIGMAcel(ijk)
!		call flush(6)
!	
!	Atualizacao do passo temporal	
!
!		3 iteracoes de Newton para raiz de SIGMAcel
!
		RSIG = SIGMAcel(ijk)
		RSIG = RSIG - (RSIG*RSIG - SIGMAcel(ijk))/(2.*RSIG)
		RSIG = RSIG - (RSIG*RSIG - SIGMAcel(ijk))/(2.*RSIG)
		RSIG = RSIG - (RSIG*RSIG - SIGMAcel(ijk))/(2.*RSIG)

!		RSIG =sqrt(SIGMAcel(ijk))
						
		DTcel(ijk) = DTcel(ijk)*(1./RSIG + .005) 				! eq. (3.38)
!		print *, "celula",ijk,"DT corretor",DTcel(ijk);call flush(6)
	enddo												! bloco ijk
!
!	Elege o novo DT para o bloco
!
	SIGMA = EPS2*maxval(SIGMAcel)				! recupera a variacao
!	SIGMA = EPS2*maximo(SIGMAcel,maxblock_size)				! recupera a variacao
	!verificar este EPS2								! relativa de Y
!     RSIG = sqrt(SIGMA/EPS2)
!	FC = 1./RSIG +.005
!print *, "SIGMA ", SIGMA
!call flush(6)

	if (SIGMA.gt.EPS.and.ITER_C.le.ITMAX) then				! Crit Convergencia do Corretor
!	 
!	Itera corretor
!		
		do ijk = 1, maxblock_size
			do j=1,nspecies
				Yp(ijk,j) = Y(ijk,j)					! Yp <- Yc
			end do
		enddo
		ITER_C = ITER_C + 1
!		print *, "iterou corretor",ITER_C
	else
!
!	Atualizacao
!												
!		DT = min(minimo(DTcel,maxblock_size),(DLTT-T))			! limita o passo temporal para nao
		DT = min(minval(DTcel),(DLTT-T))			! limita o passo temporal para nao
														! ultrapassar o limite de integracao
		T = T + DT	
print *, CONT,"T =",T,DT			
										
		ITER = ITER + 1
		ITER_C = 0
		do ijk = 1, maxblock_size
			do j=1,nspecies
				Y0(ijk,j) = Y(ijk,j)				! atualiza a composicao em todo o bloco
				P0(ijk,j) = P(ijk,j)
				G0(ijk,j) = Q(ijk,j) -  P(ijk,j)*Y(ijk,j) 
			end do
		enddo
	endif

!	CONT = CONT + 1
!IF(CONT.GE.CONTMAX) THEN STOP
end do										! while T <= DLTT
!			print *,"************************************************************************************"
!			print *, " 					Composicao em T=",DLTT
!		do ijk = 1, maxblock_size
!			print *,"************************************************************************************"
!			print *, " "
!			print *, "bloco ", ijk
!			print *, " "
!			print *, "***********************************************************************************"
!			do j=1,nspecies
!				print*,"especie",j,"Y=", Y(ijk,j)				
!			enddo
!		enddo
!stop
!print *, "saindo de aqssa"
DEALLOCATE (A)
DEALLOCATE (E)
DEALLOCATE (G0)
DEALLOCATE (P0)
DEALLOCATE (P)
DEALLOCATE (Q)
DEALLOCATE (Q0)
DEALLOCATE (R)
DEALLOCATE (Y0)
DEALLOCATE (Yp)

end SUBROUTINE aqss

SUBROUTINE cinetica(RK,Y,W,DW,P,Q,nspecies,ijkbegin,ijkend,maxblock_size,nr)
!
!	Declaracao de variaveis
!
implicit none
INTEGER, INTENT(IN) :: nspecies,nr,ijkbegin,ijkend,maxblock_size
DOUBLE PRECISION, INTENT(IN) :: RK(:,:),Y(:,:)
DOUBLE PRECISION, INTENT(INOUT) :: DW(:,:,:),P(:,:),Q(:,:),W(:,:)

	CALL rates(RK,Y,W,nspecies,ijkbegin,ijkend,maxblock_size,nr)  		        !(obtem w)
	CALL dratedc(RK,Y,DW,nspecies,ijkbegin,ijkend,maxblock_size,nr)			!(obtem dw)
	CALL fexloss(DW,P,nspecies,ijkbegin,ijkend,maxblock_size,nr)			!(obtem Pp)
	CALL fexprod(W,Q,nspecies,ijkbegin,ijkend,maxblock_size,nr)			!(obtem Qp)
	
end SUBROUTINE cinetica

double precision function Pade(x)   ! (eq. 3.29)
	double precision  x
	Pade = (180.*x*x*x + 60.*x*x + 11.*x + 1.)/(360.*x*x*x + 60.*x*x + 12.*x + 1.) 
end function 

double precision function minimo(x,n) 
!******************************************************************************
!*				   Minimo de um n-vetor nao negativo	      *
!******************************************************************************
	double precision :: x(n)
	integer :: i,n
	minimo = x(1)
     do i = 2,n
		if(x(i).lt.minimo) minimo = x(i)
	end do
end function

double precision function maximo(x,n)
!******************************************************************************
!*				   Maximo de um n-vetor nao negativo				   *
!******************************************************************************   
	double precision :: x(n)
	integer :: i,n
	maximo = x(1)
     do i = 2,n
		if(x(i).gt.maximo) maximo = x(i)
	end do
end function 

END MODULE mod_chem_spack_qssa
